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ABSTRACT 



The iGoldreich and Ward! ( 1l973l ) (axisymmetric) gravitational instability of a 



razor thin particle layer occurs when the Too mre par a meter Qt = Cpilo/yrGEp < 



1 (cp being the particle dispersion velocity). IWardI ( 119761 . 120001 ) extended this 



analysis by adding the effect of gas drag upon particles and found that even when 
Qt > 1, sufficiently long waves were always unstable. lYoudinl (l2005al Jbl) carried 
out a detailed analysis and showed that the instability allows chondrule-sized (~ 1 
mm) particles to undergo radial clumping with reasonable growth times even in 
the presence of a moderate amount of turbulent stirring. The analysis of Youdin 
includes the role of turbulence in setting the thickness of the dust layer and in 
creating a turbulent particle pressure in the momentum equation. However, he 
ignores the effect of turbulent mass diffusivity on the disturbance wave. Here we 
show that including this effect reduces the growth-rate significantly, by an amount 
that depends on the level of turbulence, and reduces the maximum intensity of 
turbulence the instability can withstand by 1 to 3 orders of magnitude. The 
instability is viable only when turbulence is extremely weak and the solid to gas 
surface density of the particle layer is considerably enhanced over minimum-mass- 
nebula values. A simple mechanistic explanation of the instability shows how the 
azimuthal component of drag promotes instability while the radial component 
hinders it. A gravito-diffusive overstability is also possible but never realized in 
the nebula models. 

Subject headings: Protoplanetary Disks; Stars: Planetary Systems 



1. Introduction 



1.1. Preliminary Remarks 



It is believed that clumping of solid material to form the terrestrial planets and the 
putative cores of gas giants involved three stages. The first and third stages are relatively 
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well understood. First, grains that survived shocked entry into the solar nebula and those 
that condensed from the cooling gas, collided and stuck by van der Waals attraction; such 
a process is coiisiderably speeded up by turbulence but i s still effective in laminar disks 
( lWeidenschilling||l980l . Il984j : iDuUemond and Dominikll2005l ) and is able to form cm-or-larger 
sized particles. Thereafter, growth by binary accretion becomes problematic in a turbulent 
nebula due to shattering. In the third stage, planetesimals, bodies 1 km and larger which are 
akin to the present day asteroids, merged in binary fashion by physical and/or gravitationa l 
capture as they collided (IKokubo and Idall2000l : IChambersll200ll : iKenyon and Bromleylll993l ) . 



Lea st understood is the middle stage, namely, the gr owth from cm to km-sized bodies. 
Early on ISafronovl (Il972l ) and lGoldreich and Ward! (119731 ) (henceforth collectively referred to 
as SGW) independently suggested that particles settled to the midplane and underwent a 
gravitational instability (GI). This caused clumps of some charac teristic size to contra c t unti l 
centrifugal force became strong enough to balance self-gravity. iGoldreich and WardI (119731 ) 
suggested that bodies of size ^0.5 km having the density of solid material could form in 
this way on a dynamical timescale. 

The SGW scenario ignored the presence of global turbulence in the solar nebula which 
stabilizes the dust layer against gravitational instability by making the basic state particle 
layer more diffuse and introducing a turbulent pressure into the particle dynamics. A pos- 
sible source of disk turbulence is magneto-rotational instability (MRI). It has been argued 
( IGammid Il996l ) that MRI may be confined to only the upper layers of a disk which are 
sufficiently ionized by cosmic rays. Nevertheless, this does not imply a completely quiescent 
mid-pl ane since vortical eddies i n the unstable layers can induce fluctuations at the mid- 
plane. iDuUemond and DominikI ( 120041 ) found that laminar settling of dust grains could not 
account for the observed flaring of disks suggesting that turbulence is keeping small particles 
aloft throughout the disk height. 

Even if one discounts global turbulence, an obstacle to the SGW picture of midplane 
GI remains: the dust disk itself can generate turbulence. The gas component feels a small 
outward directed pressure force and therefore rotates at a slightly slower speed than the 
Keplerian speed of solid particles. This creates an Ekman-like layer whose density strat- 
ification is not strong enough to stabilize the layer. The turbulent state is such that the 
critical particle density required for the mid- plane gravitational instability is not realized 
(1Weidenschilli^ll980l . Il984j : ICuzzi et al.lll993l ). 



A variant of midplane gravitational instability involves suppression of turbulence in 
the mixed particle-gas layer by the vertical density gradient. This variant is valid only 
for very small particles, which are so well trapped to the gas that the particle-gas mixture 
behaves as a single fiuid. In this limit, ISekiyal ( 119981 ) calculated the (mean) vertical structure 
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of the turbulent dust layer by assuming that the mean profiles of density and velocity lie 
at the instability/stability boundary, i.e., characterized by a Richardson number equal to 
its critical value of 1/4. He found that if the ratio of solid to gas surface density in the 
disk is much larger than the cosmic value, the layer achieves the critical density needed 
for GI. This type of analysis produces a cusp at the mi dplane in the mean den sity profile 
as the surface densi t y of s olids in the disk is increased. lYoudin and Shul (120021 ) and later 
Youdin and ChiangI (120041 ) studied this situation further, suggesting that the cusp was not 
a mathematical artifact but indicative of the inability of the turbulence to stir up particles 
near the midplane. They propose d that the requir e d enr ichment of solids could occur as 
particles drift inward by gas drag. iGaraud and Lin] ( 120041 ) studied the Richardson number 
instability as the particle layer slowly settles. They also find (their Fig. 11) that, at the levels 
of the pressure force appropriate to a minimum mass nebula, gravitational instability sets in 
before the Richardson number instability only if the column density of particles is ~ 1/5 the 
column density of the gas; this represents a factor of ~ 30 increase over the minimum-mass- 
nebula model, in agreement with previous works. Because of the small particle limitation, 
this variant of midplane GI is p recluded by even the faintest breath of global turbulence 
( ICuzzi and Weidenschillingil2006l ). 



Alternate scenarios of growth in dense midplane layers, which do not rely on grav- 
itational instabilities , have been developed by Weidenschilling in a number of papers (see 
Weidenschillingll2000l for a review). A review of this stage is giveri bv lCuzzi and Weidenschilling 
( 20061 ) and recent work ( Cuzzi et al. 2001 : Johansen et al. 2006 ) suggests that the problem 
is more complex than envisioned in the early works. 



1.2. The gas-drag mediated instability 



The subject of the present work is a gravitationally driven instability made possible by 
gas drag. What makes the instability interesting is that, although it is relat ively slow, it is 



Wardl JlOrdboOoh in the absence 



unconditional. The instability was proposed and studied by 

of turbulence and associated particle dispersion. ISafronovl (119871 Il99ll ) also briefly considers 
this instability under the same conditions and writes down the dispersion relation in the 
limit where particle stoppirig time is short compared to the time scale of the perturbation. 
A detailed analysis (lYoudinll2005al Jbl. henceforth papers I and II) suggested that the gas-drag 
assisted instability remains viable in the presence of a moderate amount of nebula turbulence. 
This analysis included the effect of turbulence in the particle momentum equation (via an 
effective pressure) and in setting the height of the disk layer. However, a turbulent diffusion 
term must also be present in the equation for particle mass conservation. Here we correct 
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this omission and find significant reductions in growth-rate. In particular, we find that the 
instabihty is not viable for nominal values of nebula turbulence unless favorable assumptions 
are made regarding local conditions. These include some combination of elevated solid/gas 
ratio, and particle size neither too small nor too large. 



During review of the revised version of this paper, we were sent a new paper, lYoudin 



(120111 ). that also corrects for the lack of turbulent mass diffusivity in the 2005 papers. In 
addition, it uses updated models for the radial turbulent diffusivity of particles, the radial 
particle dis persion, and the height of th e particle sub-disk. These models are based on the 
analysis of lYoudin and Lithwickl (j2007| ) and the corrections they imply become important 



for particles whose stopping times 4 non-dimensionalized by the orbital frequency Qq are 
such that Ts = Qts > 1. The present calculations were therefore redone to incorporate the 
newer models, however, the basic conclusions remained the same. 



2. Analysis 

2.1. Dispersion Relation 

When a mass diffusion term with coefficient D is included, the linearized and vertically 
integrated mass and momentum equations for the particles, equations (6) and (7) in paper 
I, become 

da du ^d'^a 



.2 



dt dx dx 



dv 



+ {2-q)u^^ = (3) 
at ts 

where q = 3/2. The equations are written relative to a box at radius R revolving at the local 
angular velocity RQ{R) of particles; we have defined Qq = ^{R)- Affixed to the box are 
local Cartesian coordinates {x,y), where x is radial and y is azimuthal. The corresponding 
velocity components are {u, v): they represent turbulent mean quantities that have then been 
vertically averaged. The quantity a is the relative perturbation in particle surface density 
defined so that: 

S = Sp (1 + a) , (4) 

where Sp and S are the basic state and total (basic state + perturbation) surface densities, 
respectively. Throughout the analysis used to obtain ([I])-®, one assumes that Sp is locally 
uniform and can therefore be freely moved into and out of x and y derivatives. The basic 
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state velocity of particles has been taken to be Keplerian in deriving ([I])-(l3]). The gas flow has 
also been taken to be Keplerian and is left unperturbed. In this model therefore, there is only 
one-way coupling from the gas to the particles. (One con sequence of two-way coupling are 



streaming instabilities, e.g., lYoudin and GoodmanI ( 120051 )). Figure 3 (lower) in ICuzzi et al. 



( 119931 ) shows that the turbulent mean flow departs from being Keplerian by only about 
0.2%. The basi c state flow of the particles should also have an inward radial drift arising 
from gas drag (lAdachi et al.lll976l ): this has been left out in ([I])-([3]). However, since the 
(vertically averaged) drift velocity is locally uniform with respect to x and y, it would cause 
the instability wave to merely drift inward without affecting its local growth-rate. Effects 
that involve the fact that the drift velocity has a vertical dependence are not included in 
the analysis. We found (except for the gravito-diffusive mode that was never realized in the 
nebula calculations) that the phase and group velocity of the instability wave is zero in the 
drifting frame. This comes about simply because the wave oscillation frequency is zero in 
some neighborhood of the most amplified wavenumber. Thus, the wave drift speed is the 
same as the particle drift speed and this will limit the total amount of wave growth. This is 
accounted for in the study by comparing the e-folding growth time to a local drift time scale 



ift 



R/Ud: 



rift- 



(5) 



Complete loss of solids from the entire nebula is a weaker constraint. 



Axisymmetric disturbances are being considered. Thus, the only spatial derivatives that 
appear are those with respect to the radial coordinate x. 

The quantity x = + n' consists of the gravitational potential $ and 11' = aCp which 
arises from modeling the effect of turbulence on particle momentum by an ef fective pressure, 
Cp being the radial particle dispersion velocity ( lYoudin and Lithwickl 120071 ). The effective 
pressure term arises from the rr component of the particle Reynolds stress tensor in the 
Reynolds-averaged par ticle mo n ientum equation. It should be noted that neither the present 
treatment nor that of lYoudinl (120111 ) has incorporated all the effects of turbulence in the 
particle momentum equation. A Reynolds average of the original equation reveals that, 
in addition to an effective pressure, several other turbulent correlations ar ise that require 
closure models. Compared to previous treatments ( IWardll2000l : lYoudinll2005al Jbl) the only new 
term is the right hand side of ([1]) and represents turbulent diffusion of surface density. This 
term arises from the correlation (p'^j) in the Reynolds averaged particle mass conservation 
equation (the primes denote turbulent fluctuations and the angle brackets denote a suitable 
Reynolds average). This has been modeled using gradient diffusion as —Ddp/dxj. 

Substituting into equations (II])-(l3]) perturbations of the form 



cr, u, v)e 



i(kx—ujt) 



(6) 
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where [a, u, v) are complex constants, gives a linear homogeneous system 



A(a, u, v) 



0, 



(7) 



where A is a matrix. The solvability condition for this system, namely, that the determinant 
of A vanish, gives the dispersion relation: 



where 



{oo + i/Q [C{k) -{uj + iDk^){uj + i/Q] + inliDk^ - t;^) = 0, 



C{k) = nf, + k'd - 2TiG^pkT{khp). 



(9) 



Here is the height of the particle sub-disk and T{khp) = 1/(1 + khp) is a factor that 
approximates finite thickness effects on the potential of self-gravity in the context of a ver- 
tically integrated model. In the limit tg ^ oo (vanishing drag) and -D — t- 0, ([8]) reduces 
to the Safronov-Goldreich-Ward (SGW) form = C{k). Since turbulent mass diffusivity 
occurs only in the product Dk"^, its effect is arbitrarily small for sufficiently long waves. 
However, the most amplified wave has a finite k and therefore turbulent mass diffusivity has 
a non-negligible effect on growth-rate. 



Following lYoudiru (j2005al ) let us introduce the non-dimensional variables 
7 = —ico/flo, n = 7rGSp/c/f2Q, = i^o^s, 
and the Toomre and Roche parameters 



Qk 



h VL^ 



(10) 



111 



The additional parameter resulting in the present case is the non-dimensional diffusivity: 



D = 



(7rG'Ep)2- 



(12) 



The quantity Re(7) gives the growth-rate and the parameter Qr is present only because 
finite thickness effects on self-gravity have been retained. The dispersion relation ([8]) may 
then be written as: 



where 



F{k, Qt, Qr) = 1 - + QW- 



0, 



1 + kQyi 



(13) 
(14) 



- 7- 



As can be checked, the condition F < gives instabihty in the hmit of zero gas drag (ts — )■ oo) 
and zero mass diffusivity — ?■ 0). Equation f lT3|) expands to the following cubic equation 
for the 7: 

7^ + (2r-i + Dk^'^ 72 + + F + 2Dt^T;^^ 7 + {F - I) + Dk^ (l + r'^) = 0. (15) 

Using Descartes' rule of signs and clever reasoning, lYoudinl (l2005al ) concluded that a nec- 
essary and sufficient condition for instability (for the case of zero mass diffusivity) is that 
F{k, Qt, Qr) — 1 < 0. A similar analysis does not appear to be possible in the present case. 
However, lYoudinl (l201l[ ) has shown that, even with mass diffusivity, one can always find a 
wavenumber k such that the system is unstable. 

The slow instabilities considered are such that the growth-rates can be much smaller 
than the particle stopping rate, i.e., 7 ^ I/ts. In this limit, an explicit expression is obtained 
for the growth-rate from ( !T3|) : 

7approx = rs(l - F) - Dk^ {t^ + l) , (16) 

which clearly shows the damping effect of D. The validity of ( IT6|) can be checked a posteriori. 



2.2. Implementation 

To implement the analysis for various dis k condit i ons w e follow exactly the treatment 



of papers I and II updated with the models in lYoudinl ( 120111 ) for the effect of turbulence on 



particles. Defining w = R/AU, the nebula model employed is 



Sg = 1700/gtu-3/2 gm cm-2, (17) 
Ep = 10/pti7-3/2 gm cm-2, (18) 
lO^w-^/^cm s-\ (19) 



for the gas-disk surface density, particle-layer surface density, and sound speed, respectively. 
Apart from the factors /g and /p, this is just the minimum mass model. The calculation of 
the non-dimensional stopping time, Tg, uses the Epstein and Stokes formulas as appropriate: 

^ r 4 X 10-\a/mm)w'/'f-\ a/A^fp < 2.002; 

\ 9 X 10-5(a/mm)2(c^7/0.3)~5/^ otherwise, ^ ' 



where 

-^mfp 



L-'^'-'', (21) 
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is the mean free path. The resulting Ts depends on both particle size and radial location in 
the disk. Equation fl20l) is valid when the particle Reynolds number Re^ = lat^ujv < 1. 
This condition is satisfied in all the plots we present by adjusting the range of the abscissa 
if necessary. Here Au is the magnitude of the par t icle v elocity relative to the gas which 
depends on Ts according to equation (Al) in lYoudinl (1201 ll ). and z/ = 2.45 x lO'^tu^/^/"^ cm^ 
s~^ is the kinematic viscosity. 

T he radial compo nent of the particle dispersion velocity, Cp, is calculated using equation 
(27) in lYoudinl J201lh : 

^ (l + 2r,2 + (5/4^=^)1/2 



The radial mass diffusivity due to turbulence is written as: 



ttgCg. 



D 



Sct ' 



(22) 



(23) 



where is the tur bulent rnoment um diffusivity of the gas and Sct is the turbulent Schmidt 
number. Following lYoudinl ( 120111 ) 



Sct 



:i + r: 



2\2 



1 + Ts + 4r2 ■ 

The turbulent viscosity of the gas, z/g, is defined via an a parameter: 

Z/g = OgCg/fio, 



(24) 



(25) 



where the turbulence parameter is a measure of the local turbulence intensity of the gas 
and Cg is the sound speed of the gas. After some substitutions one obtains: 



D = a. 



1 + r, + 4r2 / QTgS 



where Qxg is the Toomre parameter of the gas disk: 

Cgfio 



(26) 



For the present nebula model we obtain: 



loVp-'(t^)- 



(27) 



(2^ 



whose large leading coefficient accounts for the sensitivity to observed in the results. 
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an 



The model of I Youdiru (1201 ll ) used to determine the height hp of the particle disk is: 

hr, = ho 

where 



27rGS 



(29) 



1 + 2/Qr. 



(30) 



Equations ( l29l) and ( 130|) lead to a quadratic equation for hp/hg which is explicitly solved. 

Finally, we discuss the constraint due to inward radial drift of the instability wave as 
motivated earlier. We have 



t. 



drift 



R/\ud- 



rift I 



(31) 



The analysis of iNakagawa et al.l (119861 ) . extended to include terms quadratic in Tg, gives for 
the equilibrium (i.e., neglecting acceleration terms) radial drift speed of the particle layer: 



Wdrift = -2?7Mk 
where uk = fioR is the Keplerian speed, 

V = 



(32) 



R dp iq^in-3 1/2 
= 1.3 X 10 w ' , 

2pgUf^ oR 



(33) 



is the pressure gradient parameter, and 



p{z)/pg{z) is the ratio of particle to gas 



density (the so-called particle loading) which in general depends on the vertical position z 
in the disk. Using vertical averages for both numerator and denominator we estimate 



S hg 
Eg h 



(34) 



2.3. Mechanistic Interpretation of the Instability 

Before we begin let us note that gas drag (identified by the parameter Xg) enters not 
only as an explicit term in the momentum equations but also through setting the turbulent 
diffusivity, effective pressure, and height of the particle layer. That is, also affects the 
parameters D, Qt, and Qr. In this sub-section we focus only on the role of the explicit gas 
drag terms in the momentum equations. 

Gravitational instability occurs when the increase of mutual gravity between approach- 
ing particles outweighs restoring forces. As is well known, an important restoring force in 
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Keplerian systems and for axisymmetric perturbations comes from the outwardly increasing 
angular momentum. Consider a perturbation mode (in radial velocity and density) such that 
particles accumulate in a toroidal region and deplete in the neighboring region. Self-gravity 
is not needed at this stage. A circle of particles at the inner edge of this region is displaced 
outward. If gas drag is absent, the circle will conserve its specific angular momentum rue 
and it will slow down more steeply than the local Keplerian speed slows with radius. It will 
therefore have a deficit of centripetal acceleration relative to the pull of the central gravi- 
tator and will experience an inward restoring force. Similarly, a circle at the outer edge of 
the toroidal region will be displaced inward and will speed up faster than the local Keple- 
rian value and experience an outward restoring force. This restoring force is what makes 
Keplerian disks stable via Rayleigh's criterion and is characterized by a H-^q term on the 
right-hand-side of the dispersion relation u"^ = C{ijj) in the drag-free (SGW) case. 

Let us now consider the same picture but in the presence of the azimuthal gas drag 
term in the momentum equation. In this case, the azimuthal speed of an outward/inward 
displacing particle circle will remain closer to the local Keplerian value and there will not 
be as much of a deficit /excess of centripetal acceleration compared to the pull of the central 
gravitator. The restoring force is thus diminished. This is the mechanism by which the 
azimuthal gas drag term promotes instability. Ultimately, what has to overcome whatever 
restoring force remains is the increase in mutual gravity among the accumulating particles. 

Note that for an outwardly displacing circle, the gas flow is a tailwind relative to the 
particles and gas drag supplies them with more angular momentum. 

In actual fact, the gas revolves slightly slower than Keplerian due to pressure support. 
Hence an outwardly/inwardly displaced circle will actually experience a slightly smaller tail- 
wind/larger headwind than envisioned above. The net result is that everything that was 
described in the preceding paragraphs should be imagined to be taking place in a reference 
frame drifting inward (in the vertically averaged model). 

Consider the limit of large drag; specifically, assume that the stopping time tg is very 
short compared to the time scale l/\u}\ of the perturbation. In this case the azimuthal 
speed of the displaced circle will always remain very close to that of the local gas, namely, 
Keplerian. In this limit, therefore, the deficit /excess of angular momentum of displaced 
circles and hence the restoring tendency is completely anuUed. One can see this explicitly 
by writing the dispersion relation (|8]) for this case: 

t;^ [C{k) -{u + iOe) i t^^] + nlDk^ - Vtl t-^ = 0. (35) 

One observes that the restoring force term {Qq in C{k)) is cancelled by the last term in fl35|l . 

So far, our discussion has only considered the role of gas drag in the azimuthal momen- 
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turn equation. We now wish to study the role of the radial drag term. To keep the focus 
on the drag terms in the momentum equation, the parameters Qt, Qr, and D will be fixed 
while Ts will be varied. For simplicity, we will set Qt = Qr which corresponds to the case 
where the disk height is estimated as hp = Cp/fio- 

To begin with let us set the mass diffusivity D = and plot the non-dimensional 
growth-rate Re(7) (maximized with respect to wavenumber and the three roots of the cubic) 
versus Tg] see Figure [TJ The first case, Qt = Qr = 0.2 (Figure [T^), is unstable in the 
limit of zero explicit gas drag {t^ — )■ oo) in the classical SGW way. Observe that in this 
case, increasing overall gas drag (lowering Tg) reduces the growth-rate (solid line). If one 
artificially retains only the azimuthal gas drag (dashed line; radial drag turned off) then the 
instability is enhanced in accordance with the discussion of the preceding paragraphs. If only 
radial drag is retained (chain dotted; azimuthal drag turned off) we see that the instability 
is further weakened compared to the case with both drag terms active. We thus conclude 
that azimuthal/radial drag promotes/hinders instability. This is easy to appreciate: radial 
gas drag slows radial compressive and rarefactive particle motions. 

The second case, Qt = Qr = 2 (Figure Wp), is neutrally stable without gas drag 
(ts — 7- oo) and with D = 0, i.e., neutrally stable in the SGW limit. It is only for this case 
that one may say "gas drag assists instability." Optimal enhancement occurs for Tg = 1. 
If only azimuthal drag is present (dashed line, radial drag suppressed) then the more drag 
the better , i.e., the optimum disappears, again in agreement with the discussion in earlier 
paragraphs. As Tg decreases, radial drag diminishes the instability, while azimuthal drag 
enhances it, leading to an optimum Tg. When only radial drag is present, the maximum 
growth-rate is zero and thus is not shown. An inspection of all three roots of the cubic 
shows that this makes physical sense. In the absence of gas drag there is a zero root and 
two purely oscillating roots. In the presence of radial drag alone the two oscillating roots 
are damped, while the zero root remains zero. 

Next, let us introduce a little diffusivity, D = 1.0. The case Qt = Qr = 0.2 (SGW 
unstable) is shown in figure [It. Diffusivity does reduce growth-rates (compare with figure 
H^) as expected, however, above a certain Tg (^ 2 for this case) azimuthal drag is stabilizing. 
This is contrary to our physical explanation and it will be suggested below that a gravito- 
diffusive instability mechanism is at play. Inspection of the individual roots showed that 
this behavior always corresponded to a pair of complex conjugate 7 roots with non-zero 
oscillation frequency (overstability, Im(7) 7^ 0). This can also be described as a pair of 
growing waves propagating radially inward and outward. The bump or local maximum in 
the solid line near Tg ~ 1 shows that gas drag assisted behavior is present to the left of the 
plot. This was not the case for D = 0. That is, mass diffusivity allows the gas drag assisted 
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instability to occur at lower Qt- If we add mass diffusivity to the case [Qt = Qr = 2) 
that is gas drag assisted without diffusivity, we observe (figure [T]i) the same qualitative 
behavior to the left of the plot as in D = case (apart from an expected decrease in peak 
growth-rate). However, to the right of the plot, the gravito-diffusive mode appears for Tg 
above ~ 200, and closer inspection reveals that for Ts > 400 the diffusive case has a larger 
growth-rate that the non-diffusive case. Notice that the gravito-diffusive mode survives in 
the gas-free limit Tg — )■ oo when there is some mechanism of mass diffusivity still active to 
keep D finite. An overstability with momentu ni diffusivity (v i scosity ) is known to exist in 



the context of planetary rings (for a review see ISchmidt et al.l (120091 )). In the present case, 
we have encountered overstability with mass diffusivity. We shall not speculate on whether 
it is relevant in other astrophysical contexts. Varying Qt in the Ts — ?■ oo limit (for the razor 
thin case, Qr = 0) shows (see figure [2]) that diffusivity allows the Toomre instability to occur 
at all Qt, albeit with decreasing growth rate as increases. It is for these reasons that we 
believe this mode (which is overstable in nature) to be gravito-diffusive in nature and invite 
the reader to explain it mechanistically. 

Since D decreases with increasing in the case of stirring by a gas, the question 
arises whether the gravito-diffusive mode can be realized for particles in a solar nebula. 
To investigate this, a search was conducted in the range 10~^ < Q^g < 10"'^, 0.1 < a < 100 
cm, and 0.1 < R < 100 AU. The most amplified mode never had Im(7) 7^ 0. 

Figures [T^ and [T]F show that when D = 100 the distinction between high and low 
Qt disappears. Both now display gas drag assisted behavior except at sufficiently large Tg 
where the gravito-diffusive mode appears. 



A different mechanistic interpretation of the instability has been given by lGoodman and Pindor 



( I2OOOI ) as follows. Consider a localized clump of positive density perturbation in the particle 
layer in the form of an axisymmetric ring. At the inner edge of the ring, particles will feel 
an extra outward gravitational force and, to maintain radial equilibrium in absence of the 
gas, will need to travel at less than Keplerian speed. In the presence of gas the particles will 
therefore be energized by gas drag and drift toward the center of the ring. At the outer edge 
of the ring, particles will feel an extra inward gravitational force and will need to travel at 
faster than Keplerian speed to maintain radial equilibrium. If gas drag is now turned on, 
they will be slowed down and drift inward. 



Compared to our explanation, that of iGoodman and Pindorl ( I2OOOI ) assumes equilibrium 



between gravitational forces and centripetal acceleration, i.e., it assumes zero restoring force 
due to the angular momentum gradient. This holds in the limit (very often true and discussed 
previously) of perturbation time scale much smaller than the particle stopping time. 



- 13 - 



Finally, it is worth emphasizing that, unlike particle concentration in pressure highs 
( iBarge and Sommerial Il995l : iHaghighipour and Boss! |2003| ) which occurs solely due to gas 
drag, the present instability requires both self-gravity and gas drag. 



3. Results 

Table [1] summarizes the five models considered in paper II and repeated here with 
the mass diffusivity term added. The reference model has the following parameter values: 
particle size, a = 1 mm, enhancement/depletion factors (relative to the minimum mass solar 
nebula) of gas, /g = 1, and of particles, /p = 1. Figure [3] is analogous to Figure 4 (top) in 
paper II and shows e-fold growth times for the reference case and various weak to moderate 
values of Og. (Here and henceforth, growth times are plotted for the fastest mode among 
all three roots and all wavenumbers.) The lower set of curves (of lighter weight) are for 
zero mass diffusivity. The upper curves (heavier weight) are for the case with diffusivity. 
One concludes that including mass diffusivity increases growth times by about two orders 
of magnitude in the asteroid belt region at 3 AU. It appears that in this region turbulence 
levels must be such that ctg < 10~^ to get growth times smaller than 10^ years, a typical 
disk lifetime, and the situation becomes worse in the terrestrial planet region. However, the 
instability could perhaps play a role in the outermost regions of the disk for weak turbulence 
levels. 

For the same turbulence levels as in Figure [3l we now consider the effect of particle size 
a at i? = 3 AU; see Figure HI Growth time (heavy lines) decreases for larger particles. Also 
shown is the characteristic drift time of the instability wave which also decreases with particle 
size but more slowly than the growth time giving rise to a possibility that tgrow < ^drift- One 
observes that such a cross over occurs (for particle radii in the range of the plot) when the 
turbulence levels is sufficiently weak. Furthermore, the cross-over occurs for smaller particles 
as the turbulence diminishes. 

Paper II also considers the following four alternative models each of which makes the 
instability stronger: (i) an increase in the surface density of particles by a factor of /p = 4, 
(ii) a depletion in gas by a factor of /g = 0.1, (iii) an increase in the particle radius to a = 1 
cm, and, finally, (iv) a model designated "all" that incorporates all of these changes. 

The "all" model, in which solids are enhanced relative to gas by a factor of 40 compared 
to the minimum mass solar nebula, produces the fastest growth times of all the models: 
these are plotted in figure [5] for larger turbulence intensities ag than in the previous figure 
up to the (global) value ctg ~ 0.01 required to account for observed disk accretion rates. The 
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increase in growth time at 3 AU from the non-diffusive to the diffusive treatment is now 
just less than an order of magnitude. Thus a sohds-enhanced nebula is much less affected 
by turbulent diffusivity. This is not surprising because in the limit of vanishing gas density, 
none of the bothersome gas-related obstacles to GI will occur. For the standard turbulence 
intensity of = 0.01 the instability is not strong enough (e-fold time barely less than a 
million years). It should be noted that unlike the reference model, the "all" model produces 
a peak in growth rate at about 7 AU. 

For the same turbulence levels as in Figure Owe now consider the effect of particle size. 
Figure [6] shows growth times (heavy lines) versus particle size for the particle enriched and 
gas depleted case (/g = 0.1, /p = 4). Again, these should be compared with the thin lines 
which show the corresponding drift time. One concludes that there exists a range of particles 
sizes for which tgrow < ^drift only if the turbulence is weaker than < 10"'^. 

Paper II defines a "viable" instability as the satisfaction of three conditions. (1) e- 
folding times reasonably smaller than the disk lifetime, in particular tgrow < 10^ yrs. (2) 
Growth times less than the drift time: tgrow < i^drift where tdrift is given in equation (16) of 
paper II. (3) Wavelength of the fastest mode reasonably shorter than the radius: Af < R/2. 
Figure [7] shows the maximum value of turbulence intensity Og that is still able to produce 
a viable gravitational instability of solids. The condition that is violated when exceeds 
the maximum value is depicted using different line types: solid, dashed, and dotted for the 
three conditions, respectively. One can see that under the favorable conditions represented 
by the "all" model, the instability becomes non-viable in the terrestrial planet region when 
the turbulence level exceeds Og = 10~^, and everywhere when Og > 10~^. 

Finally, Figure |8] shows the relative error in the maximum growth rate given by the 
approximate formula (I16p valid for small growth rates compared to stopping times. The 
value plotted is the relative error in the maximum growth rate with respect to wavenumber. 
For the two cases shown, the approximate formula gives answers correct to within 10% if 

TapproxTs ^ 10 

4. Conclusion 

Under conditions of the minimum mass solar nebula at the distance of the terrestrial 
planets and the asteroid belt, the gas-drag mediated instability remains viable for chondrule 
sized particles only for very weak turbulence levels < 10~*. Even under the optimistic 
conditions of enhancement of the ratio of particle column density to gas column density by 
a factor of 40, and 1 cm particles (both of which promote the instability) the instability 
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remains viable only if < 10~^. Things look a little more optimistic further out. At 10 
AU, under standard conditions the instability becomes too slow compared to nebula lifetime 
when ttg > 10~^ while under the most optimistic conditions, viability of the instability is 
destroyed by exceedingly large wavelength for «g > 2 x 10~^. Clearly, research should be 
undertaken to pin down turbulence levels at the disk mid-plane. 
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Name 


particle size, a 


gas ratio, /g 


particle ratio, /p 


reference 


1 mm 


1.0 


1.0 


/g = 0.1 


1 mm 


0.1 


1.0 


/p = 4 


1 mm 


1.0 


4.0 


a — 1 cm 


1 cm 


1.0 


1.0 


all 


1 cm 


0.1 


4.0 



Table 1: Models 
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Fig. 1. — Study of the effect of the exphcit drag terms in the momentum equation. Growth 

rate Re(7) versus Tg at fixed Qt, Qt, and D. , both azimuthal and radial drag terms 

active; , only azimuthal drag active; , only radial drag active. 



-23- 




Fig. 2. — Maximum growth rate Re(7) versus the Toomre parameter in the drag free hmit 

with finite mass diffusivity D. The razor thin disk, = 0, is being considered. , 

5 = 0; , D = 0.1; ,5 = 1; , D = 10. 
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R[AU] 

Fig. 3. — Growth times (mimimized with respect to wavelength) for the reference model. 
Two curves are shown for each hne type. The lower set of curves (hght type) is for zero 
mass diffusivity while the upper set (heavy type) is for non-zero diffusivity. Each line type 

corresponds to a different value of as follows: , = 10~^; , Og = 10~^; , 

ccg = 10"'''. Note: Growth times above 10^ yr are not displayed. 



'^0.1 1 10 100 

a [cm] 

Fig. 4. — Effect of particle radius a on growth time (heavy lines) which should be compared 
with the corresponding characteristic drift time (thinner lines) of the instability wave. Mass 
diffusivity has been included. Minimum mass nebula (/g = 1, /p = 1) at 3 AU. The same 

turbulence strengths and line types as Figure [3] are used, namely: , Og = 10~^; , 

ag = 10-6; , ag = 10"^. 
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R[AU] 

Fig. 5. — Growth times for model "all" and different turbulence strengths ctg as follows: 

, dig = 0.01; , dig = 10~^; , ag = 10~^; , ctg = 10"^. The curves with 

heavy lines are for non-zero mass diffusivity while the thin lines are for zero mass diffusivity. 
The abscissa starts at it! = 0.2 to ensure that the particle Reynolds number < 1. 
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Fig. 6. — Effect of particle size a on growth time (heavy lines) which should be compared 
with the corresponding drift time (thinner lines). Particle enriched and gas depleted nebula 
(/g = 0.1, /p = 4) at 3 AU. The same turbulence strengths and line types as Figure [5] are 
used, namely: , a„ = 0.01; , a„ = 10 
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Fig. 7. — The maximum turbulence level ctg the instabihty can tolerate. There are two 
curves for each model. The upper curve (light weight) is for zero mass diffusivity while the 
lower one (heavy weight) is for non-zero diffusivity. The line-type denotes the condition that 

sets the maximum ag-. : nebula lifetime; : radial drift; : wavelength. The 

abscissa starts a.t R = 0.3 AU in (d) and R = 0.2 AU in (e) to ensure a particle Reynolds 
number > 1. 
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Fig. 8. — Locus of the relative error in the growth-rate (maximized over wavelength) pre- 
dicted by the approximate formula ( fT6i) as compared with the actual value (separately max- 
imized over wavelength). Each locus was obtained by varying the radius R and plotting the 

relative error versus 7approx'n3- -D 7^ for both cases. : "all" model and = 10~^; 

, /p = 4 model and Og = 10^''. 



